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Abstract 

An algorithm to compute forces at the sea bed from a finite element solution to the 
mild slope wave equation is devised in this work. The algorithm is best considered 
as consisting of two logical parts: The first is concerned with the computation of 
the derivatives to a finite element solution, given the associated mesh; the second 
is a bi-quadratic least squares fit which serves to model the sea bed locally in the 
vicinity of a node. The force at the sea bed can be quantified in terms of either lift 
and drag, the likes of Stokes' formula or traction. While the latter quantity is the 
most desireable, the direct computation of tractions at the sea bed is controversial 
in the context of the mild slope wave equation as a result of the irrotationality 
implied by the use of potentials. This work ultimately envisages a 'Monte Carlo" 
approach using wave induced forces to elucidate presently known heavy mineral 
placer deposits and, consequently, to predict the existance of other deposits which 
remain as yet undiscovered. 

Keywords: waves; sediment; Berkhoff equation; mild slope wave equation; lift; drag; 
Stokes' formula; traction; placer deposits; heavy minerals; waves; refraction; diffraction; 
reflection; interference; standing waves; resonance 



1 Introduction 

The mild slope wave equation is a model for break water diffraction, reflection and re- 
fraction which has been used with considerable success for the quantitative prediction of 
ocean dynamics in a great variety of circumstances (see Booij for limitations). The 
model is linearised, assumes the sea bed to be locally fiat, uses potential theory and there 
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is no turbulence (see Berkhoff [U, Bettess and Zienkiewicz 0, GoNSALVES |T0| ). 
Despite this, a remarkable resemblance between the geometries of some heavy mineral 
placer deposits and those of computer-generated wave height envelopes (predicted using 
the mild slope wave equation for waves moving over fairly simple, idealised bathymetries) 
is documented in Childs and Shillington P|. Wave reflection, refraction, diffraction 
and resonance would appear to have played a major concentrating role in the formation 
of these deposits. 

An algorithm to compute forces at the sea bed from a finite element solution to the 
mild slope wave equation and the associated mesh is devised in this work. Two main 
components are fundamental to the logic of the algorithm. One is concerned with the 
computation of the derivatives to a finite element solution, given the associated mesh; 
the other is a bi-quadratic least squares fit which serves to model the sea bed locally 
in the vicinity of a node. There is a considerable advantage in developing a routine 
to compute the derivatives separate from the existing code (adapting the code to an 
alternative wave model would be one example). The computation of the wave number 
using a Newton-Raphson scheme and other components essential to the algorithm are 
also discussed. 

This work ultimately envisages a "Monte Carlo" approach using wave induced forces to 
elucidate presently known heavy mineral placer deposits and, consequently, to predict the 
existance of other deposits which remain as yet undiscovered. The intention is therefore 
to use the results in an empirical or qualitative (as opposed to quantitative) manner. 



1.1 Traction and the Boundary Layer Controversy 

The flow forces at the sea bed can be quantified in terms of either lift and drag. Stokes' 
formula or traction. While the latter is most desireable in physical terms, the direct 
computation of traction at the sea bed is controversial in the context of the mild slope 
wave equation as a result of the irrotationality implied by the use of potentials and the 
consequent lack of a thorough treatment of the boundary layer. Computing the traction 
indirectly (by using the solution to the mild slope wave equation as a boundary condition 
in a model more suited to boundary layer application eg. Childs 0, 0, |0| and 0), 
though not impossible, is computationally exhorbitant. The aforementioned controversy, 
practicality and the observed negligeable effect of the pressure gradient on the mechanical 
character of fluid motion in the vicinity of the bed (Yalin [0) suggest that velocity^ 
might be the more attractive option. Stokes' formula is probably the most conventional 



option advocated by classical texts such as Landau and Lifshitz |jT^. A comparative 
study involving all four approaches is ultimately what is required. 

The traction formulae are by far the most complicated and they incorporate all the 
elements necessary for the calculation of the other quantities mentioned. Lift, drag and 
the quantities necessary to evaluate Stokes' formula are all incidental to the traction 
calculation and it is for this reason that the traction algorithm is supplied as the central 
theme to this work. 



^to which Hft and drag are squarely proportional 
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This work is also concerned with the stabihty of fairly small, sediment grains, grains 
whose threshold is presently reached at deep to intermediate wave depths where the 
orbitals are relatively small. Scaling arguments suggest that an oscillatory flow in which 
oscillations are relatively small in comparison to the wave length is a potential flow 
to first approximation. The lateral extent of the sediment deposits of interest, taken 
in conjunction with observations that the convective term is negligeable (Yalin ||17||), 
suggests a fairly uniform boundary layer may be assumed. It may therefore be possible 
to ignore the exact physics of the boundary layer at the scale on which the sediments 
of interest occur, leaving the way open for the qualitative use of a traction calculated 
directly from the solution of the mild slope wave equation. Under these circumstances 
the tractional flow driving, what is assumed to be a relatively thin and uniform boundary 
layer is what is being considered. The modelled motion for a linear sea bed would be 
that of a number of layers of fluid slapping up and down, a kind of pumping action on 
the sea bed. 



2 Stress in Terms of a Solution to the Mild Slope 
Wave Equation 

The approximated velocity potential based on the solution to the mild slope wave equation 



IS 



HxuX2, xs, t) = Re {/'^(xi, X2)e-'"*} Z{xs, h) 



where $ is the velocity potential, Re{ } indicates the real part of a complex number, 
is the finite element solution to the mild slope wave equation, Z is a function which 
describes attenuation with depth, X3 is the vertical coordinate measured from mean water 
level, h is the depth below mean water level and is a frequency. The stress tensor is 
given by the constitutive relation 

a- = -p/ + /i(Vv + (Vz;)*), 

where, in terms of the approximation ([1|), 

\ \ dx\ dxi dxi dxi I J 



■iiut 



\ dxo dxp. dx-iaxo / 
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Forthcoming sections are devoted to the modeUing and computation of these values. 



3 The Analytic Derivatives of a Finite Element So- 
lution 



The finite element method approximates a solution to a problem in a finite dimensional 
subspace F''. Thus for e F'', 

nPoint 

1=1 

where nPoint is the total number of nodes, the q's are the degrees of freedom (the 
discrete solution) and the ■0i(a;)'s are the shape functions. The local approximation on 
each element is 

nNode 
i=l 

where nNode is the number of nodes per element, the c-^^'s are the local degrees of 
freedom, the ^|^^(a;)'s are the localised shape functions and is the element in question. 
Differentiating both sides of the above equation. 



dxk ■■■dxi 



nNode Pti^lS^^ 



The problem of calculating the derivatives of a finite element solution therefore translates 
directly into one of calculating the derivatives of the localised shape functions on each 
element. These localised shape functions are defined in terms of a basis as follows 

where the {0j(^)} is the basis defined on the master element domain, Vt. In this way the 
problem can be transferred into one in terms of the master element. 



3.1 The Two— Dimensional Case 

For a two dimensional problem 

djjfl djJfl d^jjf^ d^^f &^jJf^ 
dxi ' dx2 ' dxf ' dxl ' dxidx2 
(by equation (2)). (3) 

'^Notice that the symmetry of the stress tensor is preserved when introducing the approximation. 



Qfh Qfh Q2jh g2jh Q2fh 



dxi ' dx2 ' dxi 



dxi 



dxidx2 



nNode 



Forces at the Sea Bed 



5 



Applying the chain rule the first derivative of the basis with respect to the first variable 
is 



d(t)i 
96 



dilj^^'' dxi dipf^ dx2 
dxi d^i 8x2 <96 



dxi dx2 




56' 96 _ 

8X2 J 

The first derivative of the basis with respect to the second variable is 

90i 8'tpl^^ 8xi 8ipj^^ 8x2 
W2 ^ + 



8x1 56 



8x1 8x2 



8x2 96 
■ 84^^ 



96' 96 



The second derivative of the basis with respect to the first variable is 




92 



dipl'^^ 9xi 



96^ 




84^'' 8x2 
8x2 96 

8'^4f' dxi 8x 



2 8^^^ 8'^xi 8'^ip. 



96 / 8x28x1 96 96 
dxi 8x2 



9x1 9^ 



8x1 



+ 



8x18x2 96 96 8x2 8^1 



8'^xi 9^x2 / 9xi \ / 9x2 

w w 196 J ' 



^9xi 9X2 
'96 96 



8^\ 



9xi 



8i,l 



8X2 



8'ij, 



8x\ 



(e) 



9x9 



i2,/,(e) 



8x18x2 

The second derivative of the basis with respect to the second variable is 

da 




8x2 
.96 
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d'^'iljf^ f dxi 
dxj \ 8^2 , 

^2,/,(e) 



+ 



d^ipf^ dxi 8x2 



8il)f^ 8'^xi 8'^iIj. 



+ 



8hj)i 8xi 8x2 8ipf 8'^X2 
8x2 8^1 



8x28x1 8^1 8^2 8x1 8^2 



8X18X2 8^2 8^2 



8'^Xi 8'^X2 ( 8X1 



M2 



'8x2 
M2 



^8x1 8x2 



+ 



8x1 



8X1 



8X2 

8x1 
8x1 

8x18x2 



'8x2 
M2 



The cross derivative of tiie basis is 



8''(t)i _ _d_idi^dxi 8ipf^ 8x2 
d^id^2 8^1 I 8x1 8^2 8x2 8^2 



dxi 8x1 ^ 8'^tl). 



i2./,(e) 



8x1 8x2 ^ 8ip. 



8'^xi 



8x1 8^1 8^2 8x28x1 8^2 8^1 8x1 8^18^: 
8'^4f' 9xi 8x2 , 84'^'' 5^x2 



+ 



8'^4f^ (^^2 8x2 



8x1 



8^1 di2 



+ 



8x18x2 8^1 8^2 8x2 8^18^2 



8'^xi 8'^X2 8x1 8x1 8x2 8x2 f 8x1 8x2 8x1 8x2 ' 
9^2 Wid^' WiW [Wid^^d^Wi. 



8x1 

d_4fl 
8x2 

8x1 

8x2 

8^4'^ 
8x18x2 
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Collecting the above expressions together and re-expressing them in a vector-matrix 
form, 









dX2 


d^i 








d(j)i 






dx2 








Ot,2 








d'^X2 








dil 








d'^X2 


dil 




Oil 


Oil 






d^xi 


d'^X2 


1. \ 






d^id^2 







' dxi 
' dxi 






'dx2 
' dx2 







dxi dxi dx2 dx2 



dxi dx2 
2 dxi dx2 
'dxidx2 dxidx2 



dxi 

8X2 
|2„/,(e) 



dx\ 
dxi 
dxidx2 



(4) 



It follows from equation (^) that 

'd4"^ d4"^ 8^4"^ 8^4'^ 8^ij, 



2„/,{e) 



8xi ' 8x2 ' 8x1 



2 ' 



8x2 ' 8x18x2 



This equation is the formula by which the much desired shape function derivatives are 
calculated. Substituting it into equation (BI) 



gfh Qjh Q2fh Q2fh Q2fh 



8x1 ' 8x2 ' 8x1 ' 8x2 ' 8x18x2 



nNode 



(5) 



1=1 



This equation is the formula by which the first, second and cross derivatives of a finite 
element solution to a two-dimensional problem are calculated. 

The Matrix Entries: The matrix entries may all be formulated by taking derivatives 
of the finite element mapping. Taking the opportunity to develop a systematic notation 
for the purposes of the algorithm simultaneously, 

nNode nNode 

^i(^) = E 4>k{^){xi \node k) = E shape{k,l) * eCoord{k,i) (6) 

k=l k=l 

where (xj \node k) is the zth coordinate of node /c, as is eCoord{k,i), the (/>fc(^)'s are the 
basis, as are the shape{k, l)'s. The matrix entries are calculated according to 

8x- nNode Q(j)^ nNode 

E -^{^i\nodek) = E shape{k,j + 1) * eCoord{k,i) 



8^, 

dP'Xi 

w 

8'^Xi 
dii8i2 



1, 



nNode q2 



k=l 
nNode 



7 J '~i^9 \node k) 



y^i \node k) 



E shape{k,j + 3) * eCoord{k,i) 

k=l 
nNode 

E shape{k,6) * eCoord{k,i), 



(7) 



k=l 
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where the definition of the shape{k, JYs follows from the equations above. 

The Derivatives of the Basis: Obtaining formulae for the various derivatives of the 
basis is an elementary exercise in differentiation. The resulting formulae in the particular 
instance of the 8-noded quadrilateral basis (appendix, page ^T]) are listed in the appendix 



on page ^ A combined structure-flow chart diagram of the algorithm which computes 



the derivatives of a flnite element solution is given on page 



3.2 Some Test Examples 

Finite element approximations for a number of simple, analytic surfaces were devised by 
evaluating the self same functions at the nodes of the test mesh depicted in Figure 
A comparison of the various derivatives of the approximated surface with those of the 
analytic function itself conflrmed the algorithm to be working. 

Test 1: For the surface 

f{Xi,X2) = 1, 

the flrst, second and crossed derivatives were obtained to specifled precision (approxi- 
mately 16 signiflcant flgures) at all thirteen nodes. 



^2 



(1,3) 



(3,3) 



(1,1) 



(3,1) 



Figure 2: Test Mesh 
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Mesh Associated 
/with Finite Element/ 
Solution 



iElement = 

iElement + 1 



iElement = 1 



Shape and Shape 
Function Deriv- 
atives w.r.t. ^, d(ij 



Construct Q -d 



Invert Q 



Shape Function 
Derivatives w.r.t. x 



Derivatives of Finite 
Element Solution 

« (e) 



Reference Globally 



Finite Element 
Solution 




Flow of Variables 
and Sequence 

Flow of Variables 
Only 

Sequence Only 



No 

*^ I Write Results Out / 



Figure 1: Combined Structure-Flow Chart Diagram of an Algorithm which Computes 
the Derivatives of a Finite Element Solution Analytically. 
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Test 2: For the surface 

the first, second and crossed 
mately 16 significant figures) 

Test 3: For the surface 

the first, second and crossed 
mately 16 significant figures) 



f{Xi,X2) = Xi, 

derivatives were obtained to 
at all thirteen nodes. 



f{xi, X2) = 2:^-4x1 + 3 

derivatives were obtained to 
at all thirteen nodes. 



specified precision (approxi- 



specified precision (approxi- 



4 The Various Derivatives of Z(x^, h) 



The function which describes the attenuation with depth is 

cosh(«;(/i + X3)) 



Z{x3,h) 



cosh(K/i) 



where the X3 coordinate is measured from the mean water level, h is the depth below mean 
water level, vr is the usual mathematical constant and k is defined by the non-dimensional 
dispersion relation 

— = Tanh(K/i). 

Observing that h = h{xi,X2) and k = k,{xi,X2), the first, second and cross derivatives 



are accordingly formulated in the appendix on page At the sea bed where x^ = —h: 

1 



7 



dZ 



dxi 
dZ 



cosh.{Kh) ' 
-1 / 



dh h Ok 
+ — 



dx2 
dZ 



8X2, 

d^Z 



dxi 



X3=-h 



X3=-h 



X3=-h 



X3=-h 



cosh(fi;/i) \dxi ndxi 
-1 



dh h dn 
+ — 



cosh.{Kh) \dX2 K,dX2 
0, 



1 



cosh(fi;/i) 



[l + K- 



dh 



dxi J dxi 

1 \ dh dn h f2 
+ 



+ 2k — 1 — sinh(/€/i) — 



d'z 



dxi 



X3=-h 



cosh(K/i) 



hnj dxi dxi k \k 
dhV d^h 



+ hn 



dK 
dxi 



h d'^K 



K 



dxi 



dxo 



^ r, I h { 2k — 1 — sinh(/t/i) 
dxi * 
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dxidx2 



xz=-h 



dxidxs 
d'Z 



8x28x3 

d^z 



x-i=-h 



x-i=-h 



xz=-h 



dh dn h f2 , 



K dxl 



1 



cosh(K/i) 
/ dh dn 



(1 + .=) 



dh dh 



d^h 



+ 



dh dn 



dxi dx2 dx2 dxi 



+ ^ + h^ 



cosh(Kh) 
cosh.{K,h) 
cosh.{Kh) 



dh dn 
K— h h- 



dxi 

dh 
dxo 



+ h 



dxi 

dn 

dx9. 



dK dK h d'^K 



dxi dx2 K dxidx2 



and 



5 The Nodal Values of ti(xi^ X2) and its Various Deriva- 
tives 



Calculating the wave number, k, for a given depth is standard procedure. The dispersion 
relation ^ 

— = Tanh(/t/i) 

is conventionally solved using Newton's method. The resulting iterative scheme is, 

_ J Tanh(K*/i) - 1 



K 



K 



Tanh(K»/i) + hK^{l- Tanh^(/t^/i)) 



where the superscript i denotes the successive iteration from which a given solution was 
obtained. The initial guess usually taken is 



27r 



where Aq is deep water wave-length. 



Once this has been accomplished for each of the n nodes belonging to a given element, 
there is no reason why these nodal values shouldn't be regarded as a discrete solution in 
order to determine the derivatives. Substituting into equation 



dK^ dK^ d'^K^ d'^K^ 



d'^K'^ 



dxi ' dx2 ' dx\ ' dx\ ' dxidx2 



riNode 



E ^\nodeAQ{i)r'd^{i)- 



The derivatives of k can, alternatively, be calculated by the implicit differentiation of the 
dispersion relation. Considering [Q{^)Y^ d? {^) must be calculated at each node j, the 
former method is the more efficient. 
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6 The Sea Bed at a Node 

Because nodes do not necessarily coincide with individual points of bathymmetry mea- 
surement, and for the purposes of taking derivatives, a "sea bed" needs to be interpolated 
locally. A straightforward fit of an n degree polynomial to the n data points nearest a 
node, the use of cubic splines and a local least squares fit were all considered as possible 
ways to interpolate bathymmetry between individual points of bathymetry measurement. 

The manner in which available data was collected proved to be a deciding factor in the 
final choice. While the use of cubic sphnes is fairly established in the modelling of known 
surfaces, the problem with unknown surfaces is that slope information at the "knots" 
is required. Such information is never available in the raw bathymetry data. A further 
factor to consider is that the actual data sampling intervals range anywhere from slightly, 
to highly, irregular. One advantage of the least squares method is that a large data set can 
be taken into account, even individual data points weighted according to their proximity. 

The argument against fitting an n degree polynomial exactly to the nearest n points 
in the vicinity of a given node is that the use of a high degree polynomial will result 
in a totally fictitious model in cases where the actual surface is of "lower degree" than 
the polynomial used, alternatively, where the sampling intervals are poor. Fitting a low 
degree polynomial surface could result in the use of an unrepresentative data sample. 
The solution is therefore to fit a fairly simple, low degree polynomial surface to a larger 
data set. This can be accomplished using the least squares method. A method based on 
the least absolute value of the errors is preferable in theory, of course, but not in practice. 

Bi-quadratic and bi-cubic surfaces were experimented with using the method of least 
squares. The former was decided to be the better choice. Irregular data was found to 
allow extreme cases of the "wiggle" effect in the bi-cubic case. A bi-cubic surface also 
requires a far greater, hence locally less relevant data set and its greater degree is therefore 
not necessarily an advantage. In a real-life data comparison between actual measured 
depths, the depths predicted using cubic splines and those predicted using a local, least 
squares, bi-quadratic fit, a limited inspection suggested the least squares bi-quadratic fit 
to be superior. 



6.1 The Least Squares Fit of a Bi— Quadratic Function 

A generalised bi-quadratic equation has the form 

h{x, y) = ci + C2y + c^x + c^xy + c^y"^ + cqx^ 

or when written as the dot product of two vectors, 

K^^y) = [1, X, xy, x^] • [ci, C2, C3, C4, C5, Cq]. (8) 

A least squares fit makes, what is in one sense, an optimal choice of the constants, 
Ci,C2, ■ ■ ■ ,Cq. "In one sense", in that it minimises the summed squares of the errors at 
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the data points and not the summed absolute values of these errors. The sum of the 
squares of the errors, e, is 

n 

e = ^(ci + C2yi + c^Xi + c^Xiyi + c^y^ + cqx] - Zif 

i=l 

where the Zi are the n data points located at (xj, the points to which the bi-quadratic 
equation is to be fitted. In order to minimise e with respect to the unknown constants, 



de_ 
dci 
de_ 

dC2 

de_ 
dc3 
de_ 
dc4 
de_ 
dc5 
de_ 












J2{ci + C2yi + CsXi + CiXiyi + Cr^yJ + cqxD = Y^Zi 

i=\ i=\ 
n n 

X] Vii^i + ^^yi + C3Xi + CiXiyi + c^y'l + Cqx]) = ^ yiZi 

1=1 i=l 
n n 

Xi{ci + C2yi + CsXi + CiXiyi + c^yl + c^xl) = 

i=l i=l 
n n 

^iVii^l + C2yi + C^Xi + CAXiyi + Csy,^ + CeX- ) = Y ^iVi^i 

1=1 i=l 

n n 

Y Viici + C2yi + caXi + CiXiyi + Cf,y1 + cqx]) = Y vhi 



i=l 



i=l 



Y ^lici + C22/i + C3Xi + CiXiyi + c^yl + cqx]) = Y ^l^i 

i=l 



i=l 



Re-expressing the above system of equations in vector-matrix form, 





de 





■ 1 


Vi 


Xi 


•^iyi 


y\ 






' Ci ' 




" 1 


1 


1 " 






Vi 


vl 




Xiyf 


yf 


xjyi 




C2 




2/1 


2/2 • 


• yn 






n 


Xi 


•^iyi 


xi 


xjyi 


Xiyf 


^3 




C3 




Xi 


X2 ■ 


Xn 




Z2 


E 






























i=l 




Xiyf 


xfyi 


2 2 


Xiyf 






Ci 




xiyi 


X2y2 ■ 


Xnyn 










yf 




Xzyf 


yf 


xfyf 




C5 




yf 


yl ■ 


■ yl 












4 




1 2 


4 . 




. C6 . 






x\ 


■ x^ 





o 



Therefore 



Pc = Oz, 



where P and O take their respective definitions from the previous equation. Solving for 



c = p-^Oz. 



(9) 
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Substituting this result into equation (|), 

hix,y) = [1, y, X, xy, x^] • P^Oz 

where h{x, y) is the depth modelled locally by this least squares fitted, bi-quadratic 
equation. 

6.2 The Various Derivatives of h{x^ y) 

The corresponding derivatives of the sea bed are: 

dh 

dx 
dh 

dy 
dx"^ 



dy"^ 
dxdy 

6.3 The Sea Bed Normal 

The components of the sea bed normal are: 

dh 



[0, 0, 


1, 2/, 


0, 2x\ 


V^Oz 


[0, 1, 


0, X, 


2?/, 0] 


p-^Oz 


[0, 0, 


0, 0, 


0, 2]- 


P^Oz 


[0, 0, 


0, 0, 


2, 0]- 


P^Oz 


[0, 0, 


0, 1, 


0, 0]- 


P^Oz 



Ni{x,y) = = -[0, 0, 1, y, 0, 2x] ■ p-'Oz 

dh 

N2{x,y) = -— = -[0, 1, 0, X, 2y, 0]-P-'Oz 
dy 

,r / N 9h 

N3{x,y) = TTT = 1. 



and the unit normal, 



dh 

N{x,y) 



n{x,y) 



Nix,y) II2 



A combined structure-flow chart diagram of an algorithm to model the sea bed locally 
in the vicinity of a node by way of a least squares fitted bi-quadratic can be found on 
page M. 



6.4 Some Test Examples 



Data with which to test the algorithm was generated by evaluating a few simple, analytic 
surfaces at the required number of points. A comparison of outputted bathymetries and 
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Mesh 



Bathymmetry 



iPoint = 1 



Locate the nearest 
n Data Points 



Construct P 



Construct O 



Invert P 



iPoint = iPoint + 1 



Compute the Least 
Squares Constants 
c=P Oz 



Compute the Bathy- 
mmetry, it's Deriv- 
atives and Normal 




Flow of Variables 
and Sequence 

Flow of Variables 
Only 

Sequence Only 



Yes No 

JPoint < nPointp- / Write Results Out , 



Figure 3: Combined Structure-Flow Chart Diagram of an Algorithm Used to Model the 
Sea Bed Locally in the Vicinity of a Node. 
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sea-bed normals with those of the corresponding analytic function, from which the data 
was generated, showed the algorithm to be working. 

Tests 1: The trivial cases 

h{x, y) — c, c a constant 
were used to generate the input 

a; 1 1 2 1 3 2 

?/ 1 2 1 3 1 2 . 

z c c c c c c 

The algorithm calculated both depth and normal correct to specified precision (approxi- 
mately 16 significant figures). 

Test 2: For a topography containing the arbitrarily selected bi-quadratic 

h{x, y) = + 2t/2 + 3xy + 4:X + 5y + 6 

the input generated was 



X 


1 


1 


2 


1 


3 


2 


y 


1 


2 


1 


3 


1 


2 


z 


21 


35 


31 


53 


43 


48 



The algorithm calculated depth and normal correct to specified precision (approximately 
16 significant figures). 

Test 3: A real-life data comparison was made between actual measured depths, the 
depths predicted using a local, least squares, bi-quadratic fit and those predicted using 
cubic splines. A limited inspection suggested the least squares, bi-quadratic fit to be the 
superior choice. 

The algorithm is therefore considered to adequately perform the tasks for which it was 
designed. 



7 The Traction Acting on the Sea Bed 



The surface force per unit area, exerted by the fluid and acting on the sea bed, is given 

by 

t — crn 
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where cr is the stress tensor at the sea bed and n is the unit normal to the sea bed. In 
terms of the quantities discussed and formulated so far, 
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8x28x1 j cosh(K/i) 
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where p is the pressure, ^ is the viscosity, e and i denote the usual mathematical constants, 
Cij is a frequency, t is time, is the finite element solution to the mild slope wave equation, 
X3 is the vertical coordinate measured from mean water level, h is the depth below mean 
water level (with the exception of the superscript) and k is the wave number. The 
derivatives ^ etc. denote the jr etc. derivatives formulated in Subection |0|on 

axi ' 0x2 ox ' ay 



page (the variables x and y were used in place of xi and X2 so as to avoid confusion 
with the first and second data points, (xi, yi, zi) and (^2, 1/2, 2^2) respectively). 

A structure chart of the entire algorithm to compute tractions on the sea bed from a 
solution to the mild slope wave equation is given on page IT^. 
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Raw, Ungridded 
/ Bathymmetry Data 




Mesh Associated 
with Solution. 



Model Sea Bed at 
Nodes using Local 
Least-Squares Fit 



Solve for Wave 
Number using 
Newton-Raphson 



Analytic Derivatives 
of Finite Element 
Interpolation' 



Finite Element Sol-, 
/ution to Mild-Slope/ 
' Wave Equation 



Analytic Derivatives 
of Finite Element 
Interpolation' 



Assemble Vv 
Matrix 



Attenuate Surface 
Pressure Solution 
Down to Bed 



Traction 



' Write Traction Out / 



Figure 4: Structure Chart of an Algorithm to Compute Tractions from a Solution to the 
Mild Slope Wave Equation. 



^ see Figure |T| on page ^ for detail 
^ see Figure |^ on page |15] for detail 
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8 Conclusions 



The derivatives of a finite element solution can be successfully computed on each element 
using 



or ^ ^ _ ^ ^(e) 

dxi ' dx2 ' dxi ' ' dxidx2 



nNode 

= ^t'mr'd'i^) (10) 

Xi^) i=l 



where Q(C) and d\^) are defined in equation (Q) on page |^ cf'' is the discrete solution 
on each element and nNode is the number of nodes on each element. 

A bi-quadratic least squares fit (used to model the sea bed locally in the vicinity of a 
node) can be calculated according to 

hix,y) = [1, y, X, xy, x^] • P^Oz 

where P, O are the matrices defined on page |1^ and z is the vector of known depths 
used (sample points) in the vicinity of the node in question. The various derivatives of 
this sea bed can be calculated using 



dh 

dx 
dh 

dy 

d^h 

dx"^ 

dy"^ 



dxdy 

The components of the normal are then 



[0, 0, 1, y, 0, 2x] ■ P^Oz 

[0, 1, 0, X, 2y, 0] ■ P^Oz 

[0, 0, 0, 0, 0, 2] • P^Oz 

[0, 0, 0, 0, 2, 0] P^Oz 

[0, 0, 0, 1, 0, 0] ■ p-^Oz. 



dh 

Niix,y) = = -[0, 0, 1, y, 0, 2x]-P-'Oz, 
dh 

N2{x,y) = = -[0' 1' 0, X, 2y, 0] • p-^Oz, 

. X dh 
N3{x,y) = — = 1 



and the unit normal is 



dh 
n{x,y) 



Nix,y) 
N{x,y) \\2 



A bi-quadratic least squares fit would appear to be a superior method to model the sea 
bed locally in the vicinity of a node when compared to the more conventional approach 
which involves gridding and the use of cubic splines. 



The formula to compute the traction on the sea bed is given on page |16] (in terms of 
the derivatives of a finite element solution to the mild slope wave equation and a least 
squares fitted bi-quadratic model of the sea bed in the vicinity of each node). Lift, drag 
and Stokes' formula may all be calculated from elements incidental to it. 
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9 Appendix I 

The 8— Noded Quadrilateral Basis 

The basis used in conjunction with the 8-noded quadrilateral element is: 



l(0 = 


41: 


2(0 = 


^(e? + 6)(^2'-6) 

4 


3(£) = 


J(e? + 6)(^l + 6) 


4(0 = 




5(0 = 


-^(e?-i)(^2-6) 


6(0 = 


-^(6'+ei)(ei-i) 


7(0 = 


-^(e?-i)(?2-6) 


8(0 = 


4(^Na)(ei-i) 



The Derivatives of the 8-Noded Quadrilateral Basis 

The first derivatives of the 8-noded quadrilateral basis with respect to the first variable 
are: 

1(6 + 26-266-^2) 

-6 + 66 

^(-6 + 26-266 + 6') 
^(1-6^) 

^(6 + 26 + 266 + 6') 
-6 - 66 

^(-6 + 26 + 266 - 6') 
^(6'-i)- 



shape{l, 2 
shape{2, 2 
shape{3, 2 
shape{4, 2 
shape{5, 2 
shape{Q, 2 
shape{7, 2 
shape{8, 2 



5^ 

56 

d(j)2 

56 

5(^3 

56 



56 

505 

56 

54 

56 
56 

508 

56 
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The first derivatives of the S-noded quadrilateral basis with respect to the second variable 
are: 

1 



shape{l, 3 
shape{2, 3 
shape{3, 3 
shape{4:, 3 
shape{5, 3 
shape{6, 3 
shape{7, 3 
shape{8, 3 
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(66-6). 



The second derivatives of the 8-noded quadrilateral basis with respect to the first variable 
are: 
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da 
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da 
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The second derivatives of the 8-noded quadrilateral basis with respect to the second 
variable are: 



shape{l, 5) 



^''^^ - hi n 



shape{2, 5) = -^-^ 

5^2 
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shape{3, 5) = 




- ^(1 + ^0 
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The cross derivatives of the 8-noded quadrilateral basis are: 
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The Various Derivatives of Z{x3, h) 



Since Z(x3,/i) = ^^^i^4^4^ and - = Tanh(K/i), 

COSh(KAi) K 
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